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ABSTRACT 

yA/ms. Long-term observational data have information on the magnetic cycles of active stars and that of the Sun. The changes in the 
activity of our central star have basic effects on Earth, like variations in the global climate. Therefore understanding the nature of these 
variations is extremely important. 

Methods. The observed variations related to magnetic activity cannot be treated as stationary periodic variations, therefore methods 
like Fourier transform or different versions of periodograms give only partial information on the nature of the light variability. We 
demonstrate that time-frequency distributions provide useful tools for analysing the observations of active stars. 
Results. With test data we demonstrate that the observational noise has practically no effect on the determination in the the long-term 
changes of time-series observations of active stars. The rotational signal may modify the determined cycles, therefore it is advisable 
to remove it from the data. Wavelets are less powerful in recovering complex long-term changes than other distributions which are 
discussed. Applying our technique to the sunspot data we find a complicated, multi-scale evolution in the solar activity. 

Key words. Sun: activity - methods: data analysis 



1. Introduction 

It is well known that solar cycles cannot be represented by sin- 
] gle periodic variations. Besides the 11 year quasi-periodicity 
(Schwabe 1 843 ) of solar activity lots of different rhythms were 
discovered or at least suggested by different s tudie s (see e.g. 
; Erick et al. IT9971 Eorgacs-Dajka & Borkovits I2QQ6 ). One re- 
" markable and well- stu died p eriodicity of the solar variation was 
found by Gleissberg (11966b about 70 years ago. He described 
this variability as a "long-periodic fluctuation of the sun-spot 
numbers, the period being about seven or eight sun-spot cycles". 
I Erick et al. (11997b derived the length of the Gleissberg cycle as 
. j^lOO years from wavelet analysis of a 384 years long dataset 
■ (Hoyt et al. 11993b . The non-stationarity in the long-term varia- 
tion of active stars has already been found in active stars as well 
(Olah et al. 2007). 

Time-frequency representation of a non- stationary sig- 
nal yields information about characteristics of the dataset 
in the time-frequency plane. Then transient events or fre- 
quency/amplitude changes can be monitored easily with these 
tools. In the literature, a frequently used method is the wavelet- 
map and some other techniques based on wavelets. However, 
wavelet is only one of the several time-frequency representa- 
tions, and other methods may give better results than wavelets. 
The key problem is the balance between temporal and frequency 
resolution. The frequency-time uncertainty principle gives strict 
constraints for the time duration (AO and the bandwidth (Aoj = 
2;rA/) of a signal: A^ ^ Aco > 1/2. These measures of broad- 
ness are determined by the standard deviation of frequency and 
time, and the uncertainty principle with these quantities applies 
for any time-frequency distribution. The linear transforms (such 
as wavelet or windowed Eourier- transform) have a limit for joint 
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resolution of time and frequency for local structures: for example 
in general it is not possible to determine the frequency precisely 
at high temporal resolution. Bilinear time-frequency transforma- 
tion, called kernel method (Cohen 1995 ) can break the above 
principle locally, but its cost is paid by artifacts at other frequen- 
cies. 

While wavelets and other methods based on Eourier trans- 
forms are exclusively used to analyse stellar activity, more so- 
phisticated methods have been available for decades: the Wigner 
distribution (Wigner, [19 3 2b and its extensions by Cohen (11995b . 
These methods have also been proven to be useful in the analyses 
of variable star data of chaotic origin or non- stationary behaviour 
(see Buchler et al. 1995, KoUath & Buchler 2001). In this paper 
we present a method which can be used to study changes of stel- 
lar activity cycles in time using data with yearly gaps and with 
short-term signals (rotational modulation). 

As example we reanalyzed the available sunspot and radio 
records with modern time-frequency methods to give a complete 
picture of the beating of the Sun in the period range from a few 
years to 200 years. In the companion paper (Olah et al. 12008b 
we apply the method described in the present paper to 20 active 
stars, and discuss the time variability of their cycles. 



2. Time-frequency distributions 

A simple but still powerful method is the 'short-term Eourier 
transform' (STET hereafter). In this case the dataset is multi- 
plied by a Gaussian window centred around a given epoch ^o- If 
this windowed part of the signal is then Eourier transformed, the 
frequency spectrum of the data around the time to can be gener- 
ated. By shifting the centre of the Gaussian window in time, a 
sequence of Eourier spectra is obtained, which gives the time- 
frequency representation of the data. If the sampling of is 
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dense enough, then a two-dimensional map of the energy dis- 
tribution on the time -frequency plane can be plotted. 

We demonstrate this method on the yearly sunspot numbers 
- a dataset which is well known and analysed by several authors. 
Fig.[T]displays the STFT of the solar data. On the time-frequency 
distribution map individual windowed Fourier spectra are also 
plotted with thick lines. By changing the width of the Gaussian 
window, the balance between time and frequency resolution can 
be modified. All the locations on the map, of course, satisfy the 
uncertainty principle ^ Aoj > 1/2, and the ratio of time and 
frequency resolution is constant on the map. 
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Fig. 1. Short-term Fourier transform (bottom) of the sunspot 
number time series (top). Fourier- transforms of the subsets of 
the dataset are also plotted. The colour scale of relative varia- 
tions in the distributions is presented in Fig. 7. 

While in this demonstration we used a time centred Gaussian 
window to plot the curves of Fourier transforms, in the imple- 
mentation of STFT we use windowing in the Fourier space: 



(1) 



T(t, f) = 2 F(y) exp ( - - J exp(/27r^y)Jy, 

where F(y) denotes the Fourier transform of the signal. Note, 
that the inverse Fourier transform is calculated only for the pos- 
itive frequencies. This method has several advantages. At low 
frequencies it reduces the distortion due to the components at 
negative frequencies, moreover when cr is large, the transform 
reduces to the analytic function (see Buchler and Kollath 2001 ). 
Depending on the functional form of cr(f), Eq. 1 reduces to 
STFT or wavelet. With cr = fo/a, where /o is a constant fre- 
quency, the equivalent temporal window function of STFT is 
given by 



/z(t) = exp 



\ 2a^ I 



(2) 



Usually it is a good choice to set /o to the maximum or the cen- 
tral frequency of the studied period domain. Then a ~ 1 provides 
a naturally well balanced temporal and frequency resolution. It 
can be seen from Eq. 2 that the half width of the time centered 
windowing function is a/fo. Then, if cr depends on frequency, 
the width of the window also varies with frequency. Especially 
if cr(f) = f/a then T(t, f) is equal to wavelet, with the Morlet 
kernel, since the width of the temporal windowing function de- 
creases with 1 //. 
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Fig. 2. Wavelet (top and middle) and short-term Fourier trans- 
form (STFT) distribution (bottom) of the yearly sunspot num- 
bers. 
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Note that in the standard Morlet wavelet a is fixed to 1 . The 
introduction of a free parameter (a) which controls the balance 
between frequency and temporal resolution has several advan- 
tages as already described by Kollath & Szeidl (119931) and Soon, 
Frick&Baliunas(1999). 

As can be seen from the above definitions, wavelet with 
Morlet-kernel is very similar to STFT, however, for wavelets the 
balance between temporal and frequency resolution depends on 
the frequency. This feature is demonstrated in the top panels of 
Fig. [21 The parameter (/o) is defined in such a way that for a 
given a the wavelet and STFT have the same frequency/time 
resolution at the frequency /q. The top and middle panels of 
Fig. [2] show wavelets where the lower frequency part (top) and 
the high frequency part (middle) was matched to the STFT, re- 
spectively. Resolution in time increases with frequency, since 
the same number of cycles is covered in a shorter time-base, 
at higher frequencies. The disadvantage of this representation is 
that harmonic structures, which display the same characteristics 
in STFT, are distorted. Fig. [2] thus demonstrates that while the 
time-frequency distribution in the lower panel clearly show the 
synchronised variation of the 1 1 year cycle with its double fre- 
quency harmonic, the wavelet destroys this parallel structure. 

Note that we used a diff'erent scaling of the amplitudes at 
diff'erent frequencies. In time-frequency distributions the lower 
amplitude parts are easily washed out by the contributions from 
power of other components. Far from our topics but similar in 
mathematical tools, in speech analysis it is well known that by 
the application of a pre-emphasis filter the obscured parts of 
the frequency distribution can be recovered. A multilevel pre- 
emphasis Fourier filter is applied to visualise the whole studied 
period range equally. For the solar data, relative to the highest 
amplitude ^^11 year periodicity (0.07 - 0.14 c/j) the signal in 
diff'erent frequency ranges were amplified by the following fac- 
tors (Re): Re = 1.75 for 0.0 </ < 0.05, and Re = 3.5 above 
/ = 0.16 (/ is given in c/y). In the transition regions a smooth 
change in Re is used. With this pre-emphasis filtering there is 
a factor of 3.5 diff'erence in the amplitudes to get similar grey 
levels around 11 and 5.5 years periods. 

The optimal balance of temporal and frequency resolution 
depends on the signal itself. A question arises whether it is pos- 
sible to construct a distribution, where the resolution in both co- 
ordinates is optimal, but some other disadvantage is allowed as 
its cost. The answer is a more general group of methods for sig- 
nal processing introduced by Cohen (see e.g. Cohen 1995). The 
generalised time-frequency distribution is given by 



C(t, f) = l.J J J exp(-i^t - 2n ir f - m 
(D(^, T)s\t - T 1 2) sit -h T/2)dedTd^, 



(3) 



where ^(0 is the analysed time series and 0(^, r) is the ker- 
nel of the distribution that determines the specific properties 
of the distribution. The simplest time-frequency distribution is 
the Wigner- Ville transformation with 0(^, r) = 1 . This distri- 
bution works well only for special datasets, otherwise the re- 
sulting map is heavily contaminated by cross terms of the dif- 
ferent components. A bi-Gaussian kernel provides a distribution 
which, in the limit of the kernel width, is equivalent to STFT. 
This so-called pseudo-Wigner distribution (PWD hereafter) is 
widely used in sound processing. For comparison purposes we 
also use the Choi- Williams distribution (CWD) (see Choi and 
Williams 1 19891), which applies an exponential kernel. 

The top panel of Fig. [6] shows the PWD of the yearly sunspot 
numbers. It can be seen that the resolution in both coordinates 
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Fig. 4. False alarm probability that a structure appears in STFT 
(black) and CWD (grey) due to observational noise. The ampli- 
tude of the added Gaussian noise (from left to right): 0.0025, 
0.005, 0.01, 0.02, 0.04 and 0.08. The black rectangle represents 
the amplitudes of the visible real structures in Fig.O 



is improved compared to the result with STFT displayed in 
Fig. [2] (bottom panel), however, there are artificial features visi- 
ble on the PWD distribution which have no real meaning in the 
observed signal. The reason for these cross terms is the non- 
linearity of the method: e.g., if two sinusoidal waves with dif- 
ferent frequencies coexist in the signal then false power arises 
between the two ridges representing the cycles. 



3. Preprocessing observational data 

Stellar observations in reality have lots of deficiency, like ob- 
servational noise and imperfect sampling. Then time-frequency 
methods cannot be applied directly to the observational data, 
the analyses should be preceded by a sequence of preprocess- 
ing steps. 

In the case of active stars the situation is even more compli- 
cated: diff'erent timescales of variability coexist which interfere 
with the sampling in diff'erent ways. The rotational periods fall 
into the range of days to weeks, very close to the order of the 
sampling time. On the contrary, long term changes have a char- 
acteristic time above a year, so averaging and resampling on this 
level is needed to obtain information on these components. Thus, 
the preprocessing of the data requires extreme care. 

To avoid false averages due to the under- sampling of the 
rotational component, it is advisable to whiten all observations 
with the rotational period. Because the rotational component is 
modulated due to e.g. the diff'erential rotation, prewhitening can- 
not be performed simultaneously for the whole dataset. 

The main problem is that the light-curves are unevenly sam- 
pled, and in addition, they contain seasonal gaps. Although there 
exist methods to calculate e.g. wavelets of this kind of data di- 
rectly (like the adaptive wavelet), according to our previous ex- 
periments we prefer to interpolate the data when it is possible. 
When the sampling is denser than the dominant periodicities in 
the data, then a simple moving averaging can provide a continu- 
ous signal. Next, the averaged data can be smoothed with a cu- 
bic smoothing spline (Reinsch 1967 ). This operation provides a 
continuous function S (t) based on the discrete time series ({ti, Xi) 
i = 1, ^) by minimising the integral of 5"'(0^^ with the following 



constraint: 



1 

-y(xt-s(ti)f<(Ts 

n ^ 



(4) 
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Fig. 3. Time-frequency distribution of test data. Panel A: PWD of signal T; panel B: PWD of signal T with white noise; panel 
C: PWD of signal 'IF with averaging and spline interpolation; panel D: PWD of signal 'IF after Fourier filtering of the rotational 
period in the independent observational seasons and then averaging and spline smoothing; panel E: STFT of the same signal as in 
panel D; panel F: CWD of the same signal as in panel D. 



where cr^, the smoothing factor is a free parameter, which con- 
trols the smoothness of the spline. In theory its value should be 
chosen in the confidence interval corresponding to the left side 
of Eq.|4]thus cr^ should be of the same order as the standard de- 
viation of the observational noise. In practice, we usually check 
manually the resulting spline smoothed curves to make a fine 
tuning of the averaging and spline smoothing. As a final step 
S(t) is resampled to generate an equally spaced time series. 



4, Testing the methods 

4.1. Artificial test data sets 

To test the methods for a general signal representing the cycles 
of active stars, we created artificial datasets. Signal 'F consists 
of a periodic signal with a constant frequency of 0.002 c/d, 
a component with linearly increasing frequency (chirp) with a 
modulation between / = 0.0007 and / = 0.0011 c/d and a 
wave packet at / = 0.0003 c/d with a Gaussian amplitude 
modulation:^jc/7(-(r - 2500)^2000^)). The amplitudes of the 
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constant / = 0.002 c/d signal, the chirp and the wave packet 
are 0.02, 0.04 and 0.04 repectively. 

In addition to the components of signal T, signal 'IF has a 
higher frequency variation, demonstrating the rotational period 
of active stars, and a uniformly distributed noise with amplitude 
comparable to the rotational component. To test the preprocess- 
ing and the time-frequency methods together, the test signal 'IF 
was finalised by resampling the data with the time series of the 
observations of EI Eri (see Paper II) The data were averaged with 
a time base of 90 days, then the averages were spline smoothed 
with (Ts = 0.005 and resampled with = lOdays. 

4.2. Observational noise and time-frequency distributions 

It is well known that in case of periodic or quasi-periodic signals 
the Fourier transform has a very good noise rejection property. 
The same is true for time-frequency distributions. If the obser- 
vational noise and the signal are not correlated, then the method 
provides good result even in the case of large noise. A nice ex- 
ample of this nature of time-frequency distributions is the analy- 
sis of the semi-regular star V CVn (Buchler, Kollath & Cadmus 
120041) , in which two datasets are available for the same period of 
time. One of the records consists of visual magnitude estimates 
by amateur astronomers, while the other data are professional 
measurements by photoelectric photometry. It is clear that even 
the statistical properties of the noise sources are diff'erent in the 
two cases. Figure 2 in Buchler, Kollath & Cadmus (120041) clearly 
demonstrates that the time-frequency plots are almost fully inde- 
pendent of the choice of the dataset, so the eff'ect of observational 
noise is negligible for signals with good SNR and sampling. 

We demonstrate this noise rejection property of time- 
frequency distributions on our first test time- series. The top pan- 
els of Fig. [3l display the PWD of signal 'F (A) and the same 
signal modified with noise, with large relative temporal vari- 
ance, added to it (B). The amplitude of the noise is seen on the 
top panel of the figure. While the noise has practically no ef- 
fect on the lower frequencies, the / = 0.002 c/d component is 
slightly distorted. The results with signal 'F are presented only 
for PWD, but we note that the other time-frequency methods 
(STFT, wavelet, CWD) have the same properties. 

We have tested the sunspot observations for the sensitivity to 
noise, too. White or Gaussian noise with the same amplitude as 
the maximum sunspot number does not significantly distorts the 
time-frequency plots. 



4.3. Effect of uneven sampling and noise 

Because of the averaging and spline interpolation, the eff'ect of 
observational noise is mixed with the eff'ects of sampling, we 
performed our tests with ill- sampled data. The resulting time- 
frequency distributions obtained from the gapped signal 'IF are 
displayed in the lower panels of Fig. [3l where, on panel C the 
rotational signal is included to, while on panels D, E and F it is 
filtered out from the data. Averaging and smoothing spline inter- 
polation introduce a low-pass filtering of the signal, which can 
be compensated by the pre-emphasis filtering. In these cases we 
selected Re = 4.0 for / > 0.0015. With these settings and pro- 
cessing, all the components above about 2 years of characteristic 
time are recovered well, especially in the case when the rotation 
period is filtered out. Traces of the F = 0.002c/d (P = 1.37yrs) 
component can be found, but it is strongly distorted. The rea- 
son for this distortion is, that due to the sampling, these period- 
icities cannot be exactly recovered. Except this deficiency, the 




0.0010 



0.0005 



f [c/d] 
0.0020 



0.0015 



time [d] 



0.0010 ^ 
0.0005 ^^K^ 



2000 4000 6000 8000 



yrs 
1.37 



1.83 



2.74 



5.48 



inf. 



time [d] 



Fig. 5. Time-frequency transfer function of the processing. Top: 
STFT, bottom: PWD 



main time scales of the ill- sampled data can be recovered by the 
methods described above. This result gives a warning that, be- 
cause of the gapped nature of stellar data, activity cycles around 
1-2 years cannot be easily, if at all, recovered. 

The above conclusions are valid for all the time-frequency 
distributions we tested. We display the results for the filtered 
signal 'IF for three diff'erent transformations. The main compo- 
nents are recovered with all distributions. The cross terms are 
clearly visible in PWD and CWD (especially in CWD). Closely 
spaced frequencies, together with their cross terms tend to form 
loops (or 'bubbles'), like the ones visible in all PWD figures be- 
tween t ^ 2000 and 4000 days at low frequencies. These are 
well defined signatures of structures with two diff'erent frequen- 
cies, even short living ones, and their appearance remains mostly 
invariant against changes in the preprocessing (e.g. smoothing). 
STFT with our test data clearly displays the same double fre- 
quency structure, however, according to our experience, in some 
cases PWD is superior to STFT in this property. 

To give a quantitative test on the combined eff'ect of sam- 
pling and noise in a statistical sense, we performed Monte-Carlo 
simulation of the eff'ect of observational noise on signal 'IF. 
We calculated 10"^ realisations of the added Gaussian noise se- 
quence with diff'erent amplitude (standard deviation). Then for 
each realisation we calculated the same procedure as for the 
original test data. The time-frequency distribution of the dif- 
ference of the noisy and the noiseless smoothed data was ob- 
tained and the statistics of the amplitude of the largest feature 
(peak or ridge) was estimated. Fig. |4l shows the False Alarm 
Probability FAP(A), i.e., the probability that a structure with 
amplitude larger than A appears in the transform due to obser- 
vational noise. This test was performed with diff'erent standard 
deviations ((T„) of the added Gaussian noise, and for all the time- 
frequency distributions we use. On the figure from left to right 
the curves belong to cr„= 0.0025, 0.005, 0.01, 0.02, 0.04 and 
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0.08 mags, representing a wide range of observational noise. 
The different distributions behave similarly in this test, slight 
differences have been found only. The results with STFT (black) 
and CWD (grey) are displayed on the figure, the curves obtained 
with PWD are located between the two plotted ones. It is well 
seen, that the location of the curves are scaled linearly with the 
amplitude of the added noise. Only very high noise, above the 
displayed ones, distorts this trend, as the spline interpolation be- 
comes unstable. The amplitude of the observed features of the 
distributions of the original signal is displayed by a black rect- 
angle in Fig.m It is clearly demonstrated that even high amount 
of noise is negligible compared to the real signal in the time- 
frequency distribution. We have to note that this result depends 
on the sampling and the signal content of the observations, so 
it should be repeated for datasets with different characteristics. 
The nature of noise also influences the above results. We have 
performed our tests with a broad band (white) noise. For cor- 
related, e.g. color noise, with the same power, false structures 
appear with higher probability within the frequency range of the 
noise. Since the power density increases linearly with the inverse 
of the bandwidth of the noise (with the same total power), similar 
shift is expected in the FAP(a) curves. However, as our experi- 
ence suggests, in most applications the effect of a realistic noise 
turns out to be negligible. 

During the Monte-Carlo simulation the most probable time- 
frequency distribution due to observational noise can be calcu- 
lated. For pure Gaussian noise, the expected value of the time- 
frequency distribution < T(t,f) > is the same for all time- 
frequency pairs. The averaging and spline smoothing interpo- 
lation, however, acts as a low frequency filter, so < T(t, f) > 
depends on /. Because of the uneven sampling this filtering also 
depends on the time. Therefore, < Tit, f) > of the noise compo- 
nents can be used as a combined time-frequency transfer func- 
tion of the preprocessing of the data. Fig. [5] displays this transfer 
function for the added noise of cr^ = 0.04 for STFT and PWD. 
With PWD the most probably range of spurious components are 
narrower in time than that with STFT. This test confirms the 
above findings, that in the high frequency part of the distribu- 
tion (around 1 year) the amplitude of the signal is decreased. For 
times where the sampling is very poor, in addition, an amplifica- 
tion is possible at low frequencies. This amplification is located 
on the time-frequency plane at periods which are comparable to 
the length of the longest gaps. However, according to our ex- 
perience, with a thorough interactive control of the parameters, 
this effect can be drastically reduced. In the Monte-Carlo test of 
course it is not possible to do the same fine tuning, and as results 
the structures are visible in Fig. [51 



5. Application to solar cycles 

Long-term solar records have been analysed by several authors 
to date. A review of solar cycle evolution is presented by Usoskin 
& Mursula (2003 ), with an impressive list of references on the 
subject. The long-term behaviour of the sunspot group num- 
bers were analysed using wavelet technique by Frick et al. 
(119971) who plotted the changes of the Schwabe cycle (length 
and strength) and studied the grand minima. However, there are 
still remaining features which can be revealed by the use of time- 
frequency methods. Here we present some of these characteris- 
tics from sunspot numbers and solar radio data. 
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Fig. 6. PWD of the yearly sunspot number time series (top) and 
the same of the datasets generated from the January (middle) and 
July (bottom) monthly averages 
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Fig. 7. PWD of the monthly solar radio flux data. 



5.1. Variations Related to the Schwabe Cycle 

It is well known that the sunspot number is not the best physical 
indicator of solar activity, however it is the only direct measure 
we have for a relatively long time- span. Moreover, the sunspot 
number has an intrinsic fluctuation because of the arbitrariness 
of the number of visible spots and groups. If we assume that 
the averages of e.g. the January data are independent measures 
from that of the July data (i.e., the 'noise' for the averages of 
a given month are mostly independent from that of the month 
half a year before and after), then it is possible to construct two 
independent datasets: one for the January and one for the July 
data. The comparison of the results for the two sets gives an idea 
about the errors in the time-frequency distributions. Fig. [6] (mid- 
dle and bottom panels) exhibits the time-frequency distributions 
for periods longer than 2.5-yr for both data subsets. The PWD of 
the two subsets of the data is hardly distinguishable for periods 
longer than ^ 4 years. 

In Fig. [6] (top) we present the PWD of the whole yearly av- 
eraged sunspot data. The STFT (Fig.O bottom) shows basically 
the same structure, as the PWD, but with less visibility and reso- 
lution. Also, because the STFT does not contain any false cross 
terms it validates the higher resolution results of the PWD. 

Both the Schwabe and the Gleissberg cycles can be easily 
identified in Fig. [6l In the followings we denote the instanta- 
neous frequencies of these cycles by fs and /g, respectively. 
The modulation of the Schwabe cycle is well known. The pe- 
riod of the main variation of solar activity fluctuates between 9 
and 13 years which is well seen in our results (Fig. |2] and Fig. [6]). 
However, the previous investigations neglected to check the har- 
monics of the frequency fs . Since the solar cycle is not a sinu- 
soidal variation one should find some power at 2fs and possibly 
at 3fs . Since the harmonic frequencies are synchronised to fs , 
by the investigation of these harmonics one gets an independent 
estimate of the variation of the cycle length. The inspection of 
the history of the 2fs component confirms that around 1780 the 
cycle had a period of about 9 years, then after the minimal activ- 
ity it increased to 1 1 years. 

Till the middle of the XXth century the period decreased 
to 10 years, then it made a swing to a longer period. It is no- 
ticeable that the amplitude of the 2fs component is very small 



when the Schwabe period changes. It looks as if there were some 
locked phases, with almost constant period, when the power at 
2fs is high. Another important feature related to the nature of 
the 2fs component is that its amplitude is increased for the so- 
lar cycles with high amplitude more than expected from the ra- 
tio of the maximum levels. Since the even harmonics indicate 
the asymmetry of the form of a solar cycle (symmetric wave- 
forms have only odd harmonics) - this is another indication of 
the Waldmeier eff'ect and related phenomena (see e.g. Cameron 
and Schiissler 2008 J . 

The 'bubble' like structures in time-frequency plots indicate 
that for a finite time interval two frequencies coexist close to 
each other. The power is clearly split into two frequency parts 
after 1950, and it is best visible around 1970. The frequency 
splitting is the same at all three places (Sf = 0.025 - 0.03 c/y) 
which indicates a nonlinear connection between diff'erent peri- 
odicities. The splitting is also visible in the other time-frequency 
distributions (STFT, CWD), but it is the most prominent in PWD 
as already mentioned in Section 4.3. As an additional test we 
have calculated the time-frequency distribution from the 10.7cm 
radio flux of the Sun, which is available from 1947, and is a 
good proxy for its activity. The solar radio dataset was recorded 
at the Dominion Radio Astrophysical Observatory (DRAO) at 
Penticton, Canada (Tapping & Charrois [19941 ), three times a day, 
and are given in solar flux units or Janskys. The result is plotted 
in Fig. 121 and it is well seen that the same structure is found in 
this diagram, than that of Fig.[6]from sunspot numbers; however, 
we can not resolve the Gleissberg cycle from the radio data be- 
cause of the short time base. The possible frequencies appearing 
in the splitting of the power are indicated by horizontal dotted 
lines in Fig. [6] and Fig. [71 This test strongly confirms the physi- 
cal origin of the structure that we have uncovered. To unfold all 
the regularly spaced frequencies (/g, /g+^/, k^fs and k^fs -\-Sf, 
k=l,2) at least one more frequency is needed in addition to the 
well known ones (/g and fs , i.e., the Gleissberg and the Schwabe 
frequencies, respectively). 

5.2. Long term variations in the solar cycle 

The temporal evolution of the Gleissberg cycle can also be 
seen on the time-frequency distribution of the solar data. The 
Gleissberg cycle has been found to be variable as the Schwabe 
cycle. It has two higher amplitude occurrences: first around 1800 
(during the Dalton minimum), and then around 1950. A very in- 
teresting behaviour is the continuous decrease of the frequency 
(increase of period). While near 1750 the cycle length was about 
50-yr it lengthened to approximately 130-yr by 1950. This pe- 
riod variation explains why diff'erent works give diff'erent peri- 
ods of the Gleissberg cycle. This systematic change in the fre- 
quency does not support the conclusion of Lawrence, Cadavid 
and Ruzmaikin (1995 ) that the power at long periods is a finger- 
print of a period doubling cascade. The regular pattern of fre- 
quency splitting after 1950 gives the possibility to estimate the 
value of the Gleissberg period for the last decades. Sf = 0.03 c/y 
can be precisely determined from the splitting at fs and 2fs . 
Similarly fc + Sf = 0.031 - 0.032 dy is weU defined. From 
these two ranges (^/= 0.025-0.03 and /g + df= 0.031-0.032, so 
/g < 0.007 cly) we suggest that the period of Gleissberg cycle 
has already reached a level above 140 years. This result agrees 
with the the rate of frequency decrease from 1800 to 1950. 

We extended the Schove (1955) series which lists activ- 
ity maxima, minima and amplitudes, based mostly on aurora 
records before 1610, and using optical records afterwards, by the 
envelope of the recent sunspot data. The dates and values of the 
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Fig. 8. PWD of the long term envelope of the sunpot data. 



maxima were spline smoothed and interpolated to get an equally 
sampled envelope curve. The last part of the time-frequency dis- 
tribution of the envelope displays the same structure as the low 
frequency part of the PWD of the recent data (Fig.©. Fig. [8] ex- 
hibits the smoothed Schove series together with the variation of 
the Gleissberg cycle for 500 years. 

The Gleissberg cycle had a long period during the Maunder 
minimum. After 1700 the period jumped to a lower value 
(50 years) and started a slow increase. If we accept that the 
Gleissberg cycle gives the primary indication of the long term 
behaviour of sunspot numbers, we can interpret the 5 00-y ear- 
long history of solar activity as follows. The long period 
of the Gleissberg cycle is a good indicator of an extended 
Maunder minimum, while the subsequent, sudden increase of 
the Gleissberg period (after 1700) coincides with its termina- 
tion. An increase of the amplitude of the Gleissberg cycle around 
1 800, when its period was relatively short, resulted in a shorter 
(Dalton) minimum. During the last century the Gleissberg pe- 
riod became long, but now indicating a long lasting active stage 
of the Sun. Can we derive any prediction about the possible fu- 
ture of solar cycles from the extension of this history? The an- 
swer depends on the nature of the newly found splitting of the 
frequencies. If it is a sign of the restart of the Gleissberg cycle 
at shorter period, then we can expect faster changes in the long 
term activity level, perhaps a sudden termination of the high ac- 
tivity period. If the 30 year cycle (i.e.,/G+^/ ~ 0.032, ~ 30-yr) is 
just a temporary one, then a very slow decrease of activity level 
is expected. To rule out or confirm any of the two possibilities 
we have to wait for further observations or construct a physical 
model which is able to produce all the observed behaviour and 
unfold the nature and connection of all the diff'erent periodicities. 

Unfortunately, at present we are very far from the construc- 
tion of a proper physical model of the Sun, if it is at all pos- 
sible. Bushby and Tobias (12007b discussed the possibilities of 
predicting the solar cycle through mean-field models, in which 
the cycle modulations originated from stochastic or determinis- 
tic processes. They found that good prediction even of the near- 
est cycles is impossible in both cases. However, time-series anal- 
ysis techniques, which find a mapping function for the data, may 
give better results at least for a short-term prediction. 



6. Conclusions 

1 . It is demonstrated that time-frequency distributions are use- 
ful tools for analysing nonstationary time series: the short- 
term Fourier transform, the Wigner distribution and its ex- 
tensions are more useful than wavelets. 



2. Long-term changes (cycles) of active stars can be correctly 
identified with the ill-sampled observational data, after care- 
ful preprocessing. 

3. We found an extremely complicated multi-scale evolution 
in the solar activity. All the observed features in the time- 
frequency history of the Sun should provide strong con- 
straints on modelling the solar magnetism. 

Whether it is possible to predict the long scale future of solar 
activity from the similarities of the structures that have emerged 
in the last decades and those in the early history of the Sun, is 
uncertain. Its verification would be an important future project 
both in theory and data analysis, because of the probable con- 
nection between the long term changes in our climate and solar 
activity (see Velasco & Mendoza (120071) and references therein, 
for more). 
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